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Many social, biological, and technological networks consist of a small number of highly connected 
components (hubs) and a very large number of loosely connected components (low-degree nodes). It has 
been commonly recognized that such heterogeneously connected networks are extremely vulnerable to the 
failure of hubs in terms of structural robustness of complex networks. However, little is known about 
dynamical robustness, which refers to the ability of a network to maintain its dynamical activity against local 
perturbations. Here we demonstrate that, in contrast to the structural fragility, the nonlinear dynamics of 
heterogeneously connected networks can be highly vulnerable to the failure of low-degree nodes. The crucial 
role of low-degree nodes results from dynamical processes where normal (active) units compensate for the 
failure of neighboring (inactive) units at the expense of a reduction in their own activity. Our finding 
highlights the significant difference between structural and dynamical robustness in complex networks. 

Network science has witnessed many developments in the last decade because many real-world networks 
were found to have a variety of topological structures 1 " 3 . The complexity of a network structure can be 
characterized by the connectivity properties of the interaction pathways (links) between network com- 
ponents (nodes). The degree of a node is the number of its links connected to other nodes. In terms of the degree 
distribution (the probability distribution of the degrees of all the nodes in the network), complex networks can be 
classified into homogeneous and heterogeneous networks. Homogeneous networks such as random graphs 4 and 
small-world models 5 possess binomial or Poisson degree distribution, where the degrees concentrate around the 
mean degree. Heterogeneous networks such as scale-free networks 6 have a heavy-tailed degree distribution that at 
least approximately follows a power law. Data analyses have revealed that the scale-free property is found in 
social 7 ' 8 , biological 9 " 11 , technological 12,13 , and many other types of networks 1 " 3 ' 6 . 

One of the major characteristics of heterogeneous networks is that they are highly robust against random errors 
but are extremely fragile to attacks targeted at hubs 14 . Namely, the connectivity properties of a heterogeneous 
network as a whole are not so much affected by the random removal of a fraction of the nodes as compared with a 
homogeneous network; however, they are drastically altered by the preferential removal of hubs, leading to 
network fragmentation. These properties have been theoretically studied by means of percolation theory 15 " 17 . 
The structural fragility to the preferential removal of important nodes is also observed in heterogeneous network 
models incorporating dynamical flows of physical quantities on the network 1819 . Most of these studies are 
concerned with the robustness of the network structure: the measure of network function is given by the size 
of the giant component (the largest connected subnetwork), and failure nodes are removed to induce a topological 
change. However, less attention has been paid to changes in network dynamics caused by local errors in complex 
networks of dynamical units 20 , which are particularly relevant to biological robustness 21,22 . Because most real- 
world complex networks are composed of elements with internal dynamics, it is essential to consider the 
robustness with respect to nonlinear dynamics. 

Here we focus on the dynamical robustness of complex networks, which is defined as the ability of a network to 
maintain its dynamical activity when a fraction of the dynamical components are deteriorated or functionally 
depressed but not removed. Our aim is to demonstrate that the key nodes impacting the dynamical robustness of 
heterogeneous networks are low- degree nodes. This is in strong contrast to the structural robustness which is 
largely influenced by hubs. As an example of networks consisting of dynamical units, we introduce coupled 
oscillator networks 23 " 25 which have often been used to study a variety of biological phenomena including circadian 
rhythms, synchronized neuronal firing, and spatiotemporal activity in the heart and the brain. In the mammalian 
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suprachiasmatic nucleus (SCN) producing a biological clock that 
drives circadian rhythms, normal function is sustained by networks 
of SCN neuronal oscillators 26 . The self-sustained oscillation of each 
individual SCN neuron can be damped because of some envir- 
onmental conditions 27 or age-related deterioration 28 . Coupled oscil- 
lator models are useful for exploring how circadian rhythms are 
robustly generated in a population of self- sustained and damped 
circadian oscillators 29 . In the pancreatic islets, synchronized oscillat- 
ory activity of electrically coupled beta cells is involved in pulsatile 
insulin secretion that regulates blood glucose levels. A major issue 
that has been studied using coupled oscillator models is how rhyth- 
mic insulin secretion can be robustly achieved in populations of 
heterogeneous cells that are active or silent 30 . There are many other 
dynamical cellular and multicellular processes in which oscillation 
plays important roles 31 ' 32 . Moreover, coupled oscillator models are 
also relevant to electric power networks in which network compo- 
nents like power sources should be synchronized at the same fre- 
quency 33 . It is essential to investigate how the collective dynamics of 
power networks is tolerant to local disturbances of power sources 
that may lead to a cascading blackout 34,35 . 

We consider networks consisting of N oscillator nodes coupled by 
diffusive connections with fixed strength K (see Methods section). A 
normal (active) oscillator exhibits periodic oscillation when it is 
isolated, whereas a deteriorated (inactive) oscillator settles down in 
a quiescent state after transient damping oscillation (Fig. la) 36 ' 37 . 
Unless any local deterioration occurs, the network of all active oscil- 
lators produces completely synchronized oscillation. When a pro- 
portion of the oscillators are randomly inactivated with ratio p in a 
heterogeneous network (Fig. lb), we observe phase synchronized 
oscillation with different amplitudes (Fig. lc). The individual inact- 
ive oscillator, which cannot oscillate by itself, manages to continue 
weak oscillation in the network because of continuous inputs from 
neighboring active oscillators through diffusive coupling. Instead, 
the active oscillators decrease their oscillation levels as compared 
with the isolated ones in order to achieve a balance between their 
states and those of neighboring inactive oscillators. 

Results 

Behavior in the oscillator network model. The oscillation ampli- 
tude of a node is largely dependent on the degree of the node in 
heterogeneous networks (Fig. 2a), as compared with homogeneous 
networks (Fig. 2c). When a part of the oscillators are randomly 
inactivated independently of their degrees, the expected number of 
inactive oscillators connected to a node is proportional to the degree 
of the node. Therefore, within the subpopulation of active oscillators, 
the nodes with lower degree are likely to generate larger oscillation 
amplitudes because they are less affected by the neighboring inactive 
oscillators. To measure the network activity, we introduce the order 
parameter \Z\ that represents the average of the oscillation ampli- 
tudes over all the oscillators in the phase synchronization state (see 
Methods section). As the ratio p increases, the order parameter 
decreases (Figs. 2b and 2d). Whenp exceeds a critical value p c , the 
order parameter vanishes with a loss of global oscillation. For a 
sufficiently large coupling strength K, the order parameter displays 
a second-order phase transition at p = p c 36 ' 37 . The critical ratio p c , 
namely the maximum inactivation ratio for which the network 
activity is sustainable, is regarded as a measure of the dynamical 
robustness. A higher critical value implies a more dynamically 
robust network. 

Random inactivation. In the case where oscillators are randomly 
inactivated, the critical ratio p c is analytically obtained for both 
heterogeneous and homogeneous networks (see Methods section 
and Supplementary Information). The degree of the jth oscillator 
is denoted by kj (j = 1, N) and the mean degree, by (k) = 
^2jkj/N. We assume that the behavior of a heterogeneous 
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Figure 1 | An oscillator network model, (a) The behavior of the isolated 
active (red curve) and inactive (blue curve) oscillators. The active oscillator 
exhibits periodic oscillation, whereas the inactive oscillator becomes 
quiescent after transient damping oscillation, (b) A heterogeneous network 
composed of active (red nodes) and inactive (blue nodes) oscillators. Here, 
the oscillators are randomly inactivated with ratio p = 0.3. The link density 
is d ~ 0.08 (N = 100 and (k) = 8). (c) The behavior of (1 - p)Nactive (red 
curves) and pN inactive oscillators (blue curves) in the network. The 
coupling strength is set at K = 30. The solid black curve in each panel 
indicates the mean field for each subpopulation. The inactive oscillators, 
which are not able to oscillate when isolated, manage to sustain the 
oscillatory behavior because of the diffusive interactions with neighboring 
active oscillators. All the oscillators show phase synchronization. 

network with large N is governed by two mean fields cor- 
responding to the subpopulations of active and inactive oscillators. 
In order to apply the degree-weighted mean field approximation 38 ' 39 
or the annealed network approximation 40 to each subnetwork, 
we assume that the oscillators with the same degree in the same 
subpopulation are identical. Under this assumption, the critical 
ratio p c for heterogeneous networks is derived as follows: 
het _ F(K,a)-d 



Pc 



F(K,a)-F(K,-b) 



for K>K' 



het 



(i) 



where d = (k)/(N — 1) (~ (k)/N for large AT) is the link density, K 
is the coupling strength, and a > 0 and — b < 0 are the values of 
the intrinsic parameters for the active and inactive oscillators, 
respectively (see Methods section). The function F is dependent on 
the degree distribution normalized by the system size as well as the 
coupling strength and the intrinsic parameters (see equation (6)). 
The critical coupling strength, below which p h c et = 1, is given by 
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Figure 2 | Behavior in heterogeneous and homogeneous networks, (a), (c) The upper panel shows the degrees kj of individual oscillators in (a) 
heterogeneous and (c) homogeneous networks with d ~ 0.08 (N = 1,000 and (k) = 80), i<T = 30, and p = 0.55. The lower panel shows the oscillation 
amplitudes \zj\ of the individual oscillators, (b), (d) The order parameter IZI is plotted against the inactivation ratio p for (b) heterogeneous and (d) 
homogeneous networks with d ~ 0.08 (N = 1,000 and (k) = 80). The critical ratio p cy at which I Zl reaches 0, is different between the two types of networks 
for sufficiently large K. The error bars indicating the variance for 10 network realizations are invisibly small. 



K^ et ~a/ d m i m where d min = k min /N with the minimum degree k min = 
mm.i^j^ N kj. For homogeneous networks, we further assume that the 
degrees of all the oscillators are approximated by the mean degree for 
analytical calculations. Then, the critical ratio is described as follows: 

hom _a(Kd+b) 



Pc 



for K>K' 



horn 



(2) 



(a + b)Kd 

where K^ om = a/d< K^ et . The critical ratio is invariant if the product 
of the link density and the coupling strength is kept constant. Hence, 



a more strongly or densely connected network is less robust. The 
critical ratios for homogeneous and heterogeneous networks are 
plotted against the coupling strength K and the link density d in 
Figs. 3a and 3b. The results of the analytical formulae are in good 
agreement with the corresponding numerical results. For a wide 
range of parameter values, the critical ratio for the heterogeneous 
network is larger than that for the homogeneous network. This 
property can be further confirmed by the comparison between 
Fig. 4a for heterogeneous networks and Fig. 4b for homogeneous 
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Figure 3 | Comparison of the critical ratio with respect to coupling strength ICand link density d between heterogeneous and homogeneous networks 
for random inactivation. (a) The critical ratio p c is plotted against the coupling strength Kin networks with d ~ 0.08 (N = 3000 and (k) = 240). A globally 
oscillatory state with I Zl > 0 is observed for p < p n whereas a quiescent state with I Zl = 0 is observed for p > p c . The solid and dashed black curves indicate 
the analytically obtained results in equations (1) and (2), respectively. Blue diamonds and red triangles indicate the numerically obtained results. The error 
bars indicate the variance for 10 network realizations, (b) The critical ratio p c is plotted against the link density d in networks with N = 3000 and K = 30. 
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Figure 4 | Dependence of the critical ratio on coupling strength ICand link density d. The color indicates the value of the critical ratio p c in networks with 
N= 3000. (a) Heterogeneous networks for random inactivation. (b) Homogeneous networks for random inactivation. (c) Heterogeneous networks for 
targeted inactivation of low-degree nodes, (d) Heterogeneous networks for targeted inactivation of high-degree nodes. 



networks in the parameter space of K and d. Therefore, hetero- 
geneous networks are more tolerant to random inactivation than 
homogeneous ones. This result is consistent with that based on the 
structural robustness against random removal of nodes 1 " 314 " 17 . 

Targeted inactivation. In the case where oscillators are not ran- 
domly inactivated, the mean field approximation is not appro- 
priate. Thus, we numerically computed the critical ratios for 
preferential inactivation targeted at either high-degree or low- 
degree nodes in heterogeneous networks. First, the node with the 
highest (lowest) degree among the active oscillators was inactivated. 
Then, this process was repeated until the number of inactive 
oscillators became pN. The critical ratio for targeted inactivation is 
compared with that for random inactivation in heterogeneous 
networks, as shown in Figs. 5a and 5b. Compared with the random 
inactivation case, the targeted inactivation of low- degree nodes 
makes the networks significantly vulnerable for a wide range of 
parameter values. This consequence is opposite to the structural 
fragility against removal of hubs. Figures 4a and 4c also clearly 
show the difference between the critical ratio p c for random inac- 
tivation and that for targeted inactivation of low-degree nodes. 
As explained before, the active oscillators with lower degrees are 
likely to exhibit larger oscillation amplitudes (see Fig. 2a). There- 
fore, the inactivation of low- degree nodes more dramatically de- 
creases the order parameter and thereby results in a lower critical 
ratio compared with the random inactivation. The critical coupling 
strength indicated in Fig. 5a is obtained as K^ et ' 1 = a/d max , where 
dmax = kmaJN with the maximum degree k max = m&Xi<j< N kj. 
The critical coupling strength is considerably decreased by the 
inactivation of low-degree nodes, i.e. K^ et,t <K^ et . On the other 
hand, the inactivation of the high- degree nodes (hubs) makes 
heterogeneous networks slightly more robust as compared with the 
random inactivation. This is because active oscillators with higher 
degree contribute less to the order parameter (see Fig. 2a). In this 
case, the critical coupling strength indicated in Fig. 5a is given by 
K^ et,h = a/d min ~K^ et . The difference between the critical ratio p c 



for random inactivation and that for targeted inactivation of high- 
degree nodes is subtle, as shown in Figs. 4a and 4d. 

Discussion 

We have demonstrated that heterogeneous networks are highly 
tolerant to random inactivation but extremely fragile to targeted 
inactivation of low-degree oscillators with respect to the dynamical 
robustness. The former property is consistent with the structural 
robustness of heterogeneous networks 14 . The latter property sheds 
light on the critical difference between the structural robustness 
and the dynamical robustness. These properties also hold for 
coupled oscillator models in which the coupling strengths are 
not identical but uniformly or normally distributed (see Sup- 
plementary Information). Our finding results from the diffusive 
coupling that serves to balance the states of the two connected 
nodes. Active nodes function to restore neighboring inactive 
nodes, with a resulting weakening of their dynamical activity. 
Because low-degree active nodes are affected by only a few neigh- 
boring inactive nodes, they can maintain relatively high dynamical 
activity compared with high- degree active nodes. Therefore, tar- 
geted inactivation of low-degree nodes leads to a considerable 
reduction in the network dynamical activity, making the network 
dynamically vulnerable. 

This mechanism is expected to be applicable to a wide range of 
networked systems, because diffusion in various transport phenom- 
ena are commonly found in physical, biological, and engineering 
systems. For example, electrical synapses or gap junctions that are 
widely observed in the brain 41 can be modeled by diffusive coupling 
in neural networks 42 ; in particular, rhythmicity and synchrony 
among inferior olive neurons 43 can be quantitatively described by a 
neural network model composed of oscillatory neurons coupled 
through such electrical synapses 44 . Electric power networks that are 
composed of distributed power sources including renewable energy 
sources with inverters 45 ' 46 such as wind and solar power can also be 
modeled as a coupled oscillator network 33 interconnected by diffus- 
ive coupling with active and reactive power flows. 



SCIENTIFICREPORTS | 2 : 232 | DOI: 1 0.1 038/srep00232 



4 



www.nature.com/ scientificreports \^jf 

2S 




0 10 20 30 40 50 0 0-° 2 0-° 4 006 0.08 0.1 

K d 



Figure 5 | Comparison of the critical ratio with respect to coupling strength K and link density d between random and targeted inactivation in 
heterogeneous networks, (a) The critical ratio p c is plotted against the coupling strength Kin heterogeneous networks with d ~ 0.08 (N = 3000 and (k) = 
240), d max ~ 0.3, and d m { n ~ 0.04. Blue diamonds and the black solid curve indicate the numerically and theoretically obtained critical ratios for random 
inactivation, respectively. Green circles and pink squares indicate the critical ratios for the targeted inactivation of low-degree and high-degree nodes, 
respectively. The critical coupling strength is given by K^ etJ - 3.3 for the former case and K^ h ~ 25 for the latter case, as theoretically predicted, (b) The 
critical ratio p c is plotted against the link density d in heterogeneous networks with N = 3000 and K = 30. 



Our results imply that dynamical processes play an important role 
in understanding the dynamical robustness in complex networks 
composed of interacting dynamical units. One needs to appropri- 
ately model dynamical processes and define a measure of dynamics- 
based network function to mathematically treat the dynamical 
robustness in such realworld networks as biological cellular net- 
works 11 ' 21 , metabolic networks 9 , protein networks 10 , brain networks 47 , 
and electric power networks 48 . For designing dynamically robust net- 
works and planning recovery strategies against local errors, both the 
network structure and the dynamical processes should be considered 
in detail. 

Methods 

Model equations. The network of diffusively coupled oscillators is described as 
follows: 

K N 

z j =(a j + iQ-\z j \ 2 y j +-J2 A jk{zk-z j ) for j = 1, . . . ,N, (3) 

k=i 

where AT is the number of oscillators; Zp the complex state variable of the jth oscillator; 
ccp the intrinsic parameter of the jth oscillator; Q, the natural frequency; and K, the 
coupling strength. The single oscillator without coupling (K = 0), called the Stuart- 
Landau oscillator 24 , represents the normal form of dynamical systems that describe 
the nonlinear dynamics near the supercritical Hopf bifurcation at ccj = 0 49 . The single 
oscillator exhibits periodic (active) oscillation for a positive value of ccp whereas it 
becomes quiescent (inactive) for a negative value of a ; (see Fig. la). The set of node 
indices for the active oscillators is denoted by S A and that for the inactive oscillators, 
by Sj. We set otj = a > 0 for j <e S A and a ; = — b < 0 for j <e S/ 3637 . The parameter 
values are fixed at a = l,b = 3, and Q = 3. The adjacency matrix A = (Aj k ) represents 
the network connectivity, where A 



,j k = A k j = 1 if the jth and the kth oscillators are 
connected and Aj k = A k j = 0 otherwise. The degree of the jth oscillator is given by 
kj = J2k=i^jk- The homogeneous and heterogeneous networks are given by the 
Erdos-Renyi random graph 4 and the Barabasi- Albert scale-free model 6 , respectively. 



Order parameter and critical ratio. The macroscopic oscillation level of the entire 
network is evaluated by the order parameter \Z\ , where Z = ( 1 /N) Ylf= i z j- When all 
the oscillators stop oscillating, the order parameter vanishes. As the ratio p of inactive 
oscillators increases from zero, the order parameter decreases and vanishes at a 
critical ratio p = p c (see Figs. 2b and 2d). In simulation experiments, the model 
equation (3) with random initial conditions was numerically integrated by the fourth- 
order Runge-Kutta method with time step 0.1, and the order parameter \Z\ was 
calculated at t = 50000. We considered the critical value p c to be the value of p at 
which the value of \Z\ at t = 50000 falls below 10" 6 as p is increased. 

Mean field approximation. The critical ratio for random inactivation is analytically 
obtained by the mean field approximation 38,39 or the annealed network 
approximation 40 . The sum of contributions to the jth oscillator from the connected 
oscillators in Eq. (3) is represented by the local field h } = Ylk= i A jk z k- For sufficiently 
large N, the number of active oscillators in the neighborhood of the jth oscillator is 
expected to be (1 — p)kj and that of inactive ones, pkp Therefore, based on numerical 



observations (see Supplementary Information), we approximate the local field as 
follows: 

h j (t)^(l-p)k j H A (t)+pk j H I (t), 

where the degree-weighted mean fields H A (t) and Hj(t) for active and inactive sub- 
populations are given by 



H A (t) = — ^ and Hj(f) ~ 



Accordingly, the original equation (3) can be approximated as follows: 

Z) = (a, + iQ~ \ Zj \ 2 ) Zj + ^ ((1 -p)H A (t) +pH I {t) - Zj ) . 



(4) 



(5) 



The coupling term, which is dependent on the connectivity matrix Aj k in the original 
equation, is now only dependent on the degree of the node. Once the mean fields H A 
and Hf are provided, the steady oscillatory state of the jth oscillator is obtained from 
the reduced form (5). For the mean field approximation to be self- consistent, the 
mean fields calculated by equation (4) from the steady oscillatory states must be 
equivalent to the originally given mean fields. From a self-consistency analysis (see 
Supplementary Information), we derive the critical ratio (1) with 

1 N d 2 

F(K,a) = -y* - 3 , , (6) 

where dj = kj/N is the degree of the jth oscillator, normalized by the system size. 

For homogeneous networks, we further assume that the degrees of all the oscilla- 
tors are the same as the mean degree, i.e. kj = (k) = (l/N)Z ; /c ; - for all j. Then, because dj 
= d = (k)/N Tor all j, equation (6) is reduced to F(K, a) = d 2 /(d—a/K). By substituting 
this form into equation (1), we obtain equation (2). See Supplementary Information 
for further details. 
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